iSSF Vegetation Model Plotting
Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain vegetation categories as a covariate, with interactions for movement parameters.
Load packages
Import data
New names:
Rows: 1062908 Columns: 42
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(4): pheno_start, pheno_end, season, day_period dbl (35): ...1, id, burst_,
x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr... lgl (1): case_ dttm (2):
t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
Prepare the data
Code
# scale the covariates
ssf_dat <- ssf_dat %>% mutate(srtm_start_scaled = scale(srtm_start),
srtm_end_scaled = scale(srtm_end),
water_dist_start_scaled = scale(water_dist_start),
water_dist_end_scaled = scale(water_dist_end),
ndvi_start_scaled = scale(ndvi_start),
ndvi_end_scaled = scale(ndvi_end),
nbr_start_scaled = scale(nbr_start),
nbr_end_scaled = scale(nbr_end),
sl_scaled = scale(sl_),
log_sl_scaled = scale(log_sl_),
cos_ta_scaled = scale(cos_ta_)
)
# get scaling factors
sl_scale_sd <- attr(ssf_dat$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_dat$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_dat$cos_ta_scaled, "scaled:scale")
# rename the pheno categories to set a different reference level
ssf_dat$pheno_end <- relevel(as.factor(ssf_dat$pheno_end), ref = "shrub_grass")
ssf_dat$pheno_start <- relevel(as.factor(ssf_dat$pheno_start), ref = "shrub_grass")
unique(ssf_dat$pheno_start)[1] shrub_grass dry_grassland tree_grass bare_ground open_forest
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass
[1] shrub_grass dry_grassland tree_grass bare_ground open_forest
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass
Plot the used-available distributions
Code
# checking for all individuals
buffalo_ids <- unique(ssf_dat$id)
# inspect the data for dry season
ssfdat_dry %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season", subtitle = "All buffalo")+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Inspect the data for each buffalo
Code
for(i in 1:length(buffalo_ids)) {
print(ssfdat_dry %>% filter(id == buffalo_ids[i]) %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Dry season", subtitle = paste0("Buffalo ", buffalo_ids[i])) +
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
)
}`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
ssfdat_wet %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Wet season", subtitle = "All buffalo")+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
For each individual
Code
for(i in 1:length(buffalo_ids)) {
print(ssfdat_wet %>% filter(id == buffalo_ids[i]) %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Wet season", subtitle = paste0("Buffalo ", buffalo_ids[i]))+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
)
}`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Read in the saved model objects
Model summaries
Model fitting - interacting with movement
Code
tic()
buffalo.tmp <- glmmTMB(case_ ~
pheno_end +
sl_scaled +
log_sl_scaled +
cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id),
family = poisson(),
data = ssfdat_dry,
doFit=FALSE,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA)))
# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.dry.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)
summary(buffalo.pheno.dry.mvmt) Family: poisson ( log )
Formula:
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id)
Data: ssfdat_dry
AIC BIC logLik -2*log(L) df.resid
1161944 1162171 -580952 1161904 620303
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1e+06 1000
Number of obs: 620323, groups: step_id, 56393
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.943079 4.211029 -0.699 0.484616
pheno_endbare_ground 0.207230 0.016675 12.428 < 2e-16
pheno_enddry_grassland -0.004613 0.016601 -0.278 0.781098
pheno_endopen_forest -0.028587 0.020086 -1.423 0.154663
pheno_endtree_grass -0.082072 0.018900 -4.342 1.41e-05
sl_scaled 0.017588 0.008004 2.197 0.027998
log_sl_scaled -0.009450 0.008218 -1.150 0.250170
cos_ta_scaled -0.001376 0.006581 -0.209 0.834424
sl_scaled:pheno_startbare_ground -0.032633 0.017720 -1.842 0.065540
sl_scaled:pheno_startdry_grassland 0.054589 0.016693 3.270 0.001075
sl_scaled:pheno_startopen_forest -0.037419 0.017003 -2.201 0.027760
sl_scaled:pheno_starttree_grass 0.010700 0.018608 0.575 0.565258
log_sl_scaled:pheno_startbare_ground -0.041416 0.015747 -2.630 0.008537
log_sl_scaled:pheno_startdry_grassland -0.095950 0.017085 -5.616 1.95e-08
log_sl_scaled:pheno_startopen_forest 0.058683 0.017622 3.330 0.000868
log_sl_scaled:pheno_starttree_grass 0.006182 0.019581 0.316 0.752240
cos_ta_scaled:pheno_startbare_ground -0.046250 0.012817 -3.608 0.000308
cos_ta_scaled:pheno_startdry_grassland 0.009413 0.013829 0.681 0.496068
cos_ta_scaled:pheno_startopen_forest -0.036834 0.013465 -2.735 0.006229
cos_ta_scaled:pheno_starttree_grass -0.037956 0.014868 -2.553 0.010685
(Intercept)
pheno_endbare_ground ***
pheno_enddry_grassland
pheno_endopen_forest
pheno_endtree_grass ***
sl_scaled *
log_sl_scaled
cos_ta_scaled
sl_scaled:pheno_startbare_ground .
sl_scaled:pheno_startdry_grassland **
sl_scaled:pheno_startopen_forest *
sl_scaled:pheno_starttree_grass
log_sl_scaled:pheno_startbare_ground **
log_sl_scaled:pheno_startdry_grassland ***
log_sl_scaled:pheno_startopen_forest ***
log_sl_scaled:pheno_starttree_grass
cos_ta_scaled:pheno_startbare_ground ***
cos_ta_scaled:pheno_startdry_grassland
cos_ta_scaled:pheno_startopen_forest **
cos_ta_scaled:pheno_starttree_grass *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3.878 sec elapsed
Code
Code
tic()
buffalo.tmp <- glmmTMB(case_ ~
pheno_end +
sl_scaled +
log_sl_scaled +
cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id),
family=poisson(),
data = ssfdat_wet,
doFit=FALSE,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA)))
# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.wet.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)
summary(buffalo.pheno.wet.mvmt) Family: poisson ( log )
Formula:
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id)
Data: ssfdat_wet
AIC BIC logLik -2*log(L) df.resid
827723.7 827943.8 -413841.9 827683.7 442565
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
step_id (Intercept) 1e+06 1000
Number of obs: 442585, groups: step_id, 40235
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.078387 4.985386 -0.617 0.5369
pheno_endbare_ground 0.524979 0.018049 29.086 < 2e-16
pheno_enddry_grassland 0.122161 0.020143 6.065 1.32e-09
pheno_endopen_forest 0.188558 0.025055 7.526 5.24e-14
pheno_endtree_grass 0.056589 0.022771 2.485 0.0130
sl_scaled -0.002532 0.008825 -0.287 0.7742
log_sl_scaled 0.150044 0.011151 13.455 < 2e-16
cos_ta_scaled 0.084601 0.007900 10.709 < 2e-16
sl_scaled:pheno_startbare_ground -0.021958 0.018241 -1.204 0.2287
sl_scaled:pheno_startdry_grassland 0.006751 0.021412 0.315 0.7525
sl_scaled:pheno_startopen_forest -0.091184 0.022869 -3.987 6.68e-05
sl_scaled:pheno_starttree_grass -0.040299 0.025016 -1.611 0.1072
log_sl_scaled:pheno_startbare_ground -0.085673 0.018978 -4.514 6.35e-06
log_sl_scaled:pheno_startdry_grassland -0.050208 0.025023 -2.006 0.0448
log_sl_scaled:pheno_startopen_forest -0.123207 0.022896 -5.381 7.40e-08
log_sl_scaled:pheno_starttree_grass -0.160890 0.024680 -6.519 7.07e-11
cos_ta_scaled:pheno_startbare_ground -0.105115 0.013989 -7.514 5.73e-14
cos_ta_scaled:pheno_startdry_grassland -0.014564 0.017465 -0.834 0.4043
cos_ta_scaled:pheno_startopen_forest -0.159377 0.017228 -9.251 < 2e-16
cos_ta_scaled:pheno_starttree_grass -0.142022 0.018050 -7.868 3.60e-15
(Intercept)
pheno_endbare_ground ***
pheno_enddry_grassland ***
pheno_endopen_forest ***
pheno_endtree_grass *
sl_scaled
log_sl_scaled ***
cos_ta_scaled ***
sl_scaled:pheno_startbare_ground
sl_scaled:pheno_startdry_grassland
sl_scaled:pheno_startopen_forest ***
sl_scaled:pheno_starttree_grass
log_sl_scaled:pheno_startbare_ground ***
log_sl_scaled:pheno_startdry_grassland *
log_sl_scaled:pheno_startopen_forest ***
log_sl_scaled:pheno_starttree_grass ***
cos_ta_scaled:pheno_startbare_ground ***
cos_ta_scaled:pheno_startdry_grassland
cos_ta_scaled:pheno_startopen_forest ***
cos_ta_scaled:pheno_starttree_grass ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2.493 sec elapsed
Code
Plot both models
Create manuscript figures
Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
"Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass",
"sl", "log(sl)", "cos(ta)",
"sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
"log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
"cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
),
"Estimate" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"],
"SE" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Std. Error"]
) %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
if (coef_df$pvalue[x] <= 0.001){
return("***")
} else if (coef_df$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))
# Specify the order in which the coefficients should be plotted
order <- c(
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass",
"sl",
"log(sl)",
"cos(ta)",
"sl:Floodplain",
"sl:DryGrassland",
"sl:OpenForest",
"sl:TreeGrass",
"log(sl):Floodplain",
"log(sl):DryGrassland",
"log(sl):OpenForest",
"log(sl):TreeGrass",
"cos(ta):Floodplain",
"cos(ta):DryGrassland",
"cos(ta):OpenForest",
"cos(ta):TreeGrass"
)
# Specify the order in which the coefficients should be plotted
order_no_mvmt <- c(
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass"
)
# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))
coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Plot the coefficients
Code
min_x <- -0.4
max_x <- 0.65
# coefficients on the x-axis
dry_mvmt_plot <- ggplot(data = coef_df,
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 12.5, yend = 12.5,
colour = "gray80", lty = 1, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 15.5, yend = 15.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
# drop the filter argument when plotting all covariates
geom_errorbarh(data = coef_df2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
# ggtitle("Dry season") +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_mvmt_plotWarning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).
Drop the filter argument when plotting all covariates
Code
dry_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")),
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
# drop the filter argument when plotting all covariates
geom_errorbarh(data = coef_df2_no_mvmt,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
# scale_y_discrete(limits = rev(order)) +
scale_y_discrete(limits = rev(order_no_mvmt)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_mvmt_plotWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
"Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass",
"sl", "log(sl)", "cos(ta)",
"sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
"log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
"cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
),
"Estimate" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"],
"SE" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Std. Error"]
) %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
if (coef_df$pvalue[x] <= 0.001){
return("***")
} else if (coef_df$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))
# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))
coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Plot the coefficients
Code
# coefficients on the x-axis
wet_mvmt_plot <- ggplot(data = coef_df,
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 12.5, yend = 12.5,
colour = "gray80", lty = 1, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 15.5, yend = 15.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
geom_point(shape = 3, size = 2.5) +
geom_errorbarh(data = coef_df2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_mvmt_plotWarning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).
Drop the filter argument when plotting all covariates
Code
wet_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")),
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
geom_errorbarh(data = coef_df2_no_mvmt,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order_no_mvmt)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(-2.5, max_x, 0.5)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_mvmt_plotUpdating movement parameters, with vegetation interactions
Code
coef_df_dry_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
"scaled_coefficient" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"])
# step length coefficients
sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate))
sl_dry_mvmt_mvparams <- sl_dry_mvmt_mvparams %>% mutate(
scale_factor = sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# log step length coefficients
log_sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate))
log_sl_dry_mvmt_mvparams <- log_sl_dry_mvmt_mvparams %>% mutate(
scale_factor = log_sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# turning angle coefficients
cos_ta_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate))
cos_ta_dry_mvmt_mvparams <- cos_ta_dry_mvmt_mvparams %>% mutate(
scale_factor = cos_ta_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
tentative_mean_sl <- tentative_shape * tentative_scale
# shrub_grass
shrub_grass_mvmt_dry <- data.frame(
"pheno" = "ShrubGrass",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1],
"updated_scale" = 1/((1/tentative_scale) - sl_dry_mvmt_mvparams$nat_coef[1]),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1]
)
shrub_grass_mean_sl_dry <- shrub_grass_mvmt_dry$updated_shape * shrub_grass_mvmt_dry$updated_scale
# floodplain
floodplain_mvmt_dry <- data.frame(
"pheno" = "Floodplain",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[2],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[2])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[2]
)
floodplain_mean_sl_dry <- floodplain_mvmt_dry$updated_shape * floodplain_mvmt_dry$updated_scale
# dry_grassland
dry_grassland_mvmt_dry <- data.frame(
"pheno" = "DryGrassland",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[3],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[3])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[3]
)
dry_grassland_mean_sl_dry <- dry_grassland_mvmt_dry$updated_shape * dry_grassland_mvmt_dry$updated_scale
# open_forest
open_forest_mvmt_dry <- data.frame(
"pheno" = "OpenForest",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[4],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[4])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[4]
)
open_forest_mean_sl_dry <- open_forest_mvmt_dry$updated_shape * open_forest_mvmt_dry$updated_scale
# tree_grass
tree_grass_mvmt_dry <- data.frame(
"pheno" = "TreeGrass",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[5],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[5])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[5]
)
tree_grass_mean_sl_dry <- tree_grass_mvmt_dry$updated_shape * tree_grass_mvmt_dry$updated_scaleCreating Gamma distributions
Code
sl_dry_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))
# Tentative distribution
sl_dry_mvmt$Tentative <- dgamma(x = sl_dry_mvmt$x,
shape = tentative_shape,
scale = tentative_scale)
# Shrub grass
sl_dry_mvmt$ShrubGrass <- dgamma(x = sl_dry_mvmt$x,
shape = shrub_grass_mvmt_dry$updated_shape,
scale = shrub_grass_mvmt_dry$updated_scale)
# Bare ground
sl_dry_mvmt$Floodplain <- dgamma(x = sl_dry_mvmt$x,
shape = floodplain_mvmt_dry$updated_shape,
scale = floodplain_mvmt_dry$updated_scale)
# Dry grassland
sl_dry_mvmt$DryGrassland <- dgamma(x = sl_dry_mvmt$x,
shape = dry_grassland_mvmt_dry$updated_shape,
scale = dry_grassland_mvmt_dry$updated_scale)
# Open forest
sl_dry_mvmt$OpenForest <- dgamma(x = sl_dry_mvmt$x,
shape = open_forest_mvmt_dry$updated_shape,
scale = open_forest_mvmt_dry$updated_scale)
# Tree grass
sl_dry_mvmt$TreeGrass <- dgamma(x = sl_dry_mvmt$x,
shape = tree_grass_mvmt_dry$updated_shape,
scale = tree_grass_mvmt_dry$updated_scale)
# Pivot from wide to long data
sl_dry_mvmt <- sl_dry_mvmt %>%
pivot_longer(cols = -x)
# Specify the order in which the coefficients should be plotted
order <- c(
"ShrubGrass",
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass",
"Tentative"
)
sl_dry_mvmt$name <- factor(sl_dry_mvmt$name, levels = order)
mean_values <- c("ShrubGrass" = round(shrub_grass_mean_sl_dry, 1),
"Floodplain" = round(floodplain_mean_sl_dry, 1),
"DryGrassland" = round(dry_grassland_mean_sl_dry, 1),
"OpenForest" = round(open_forest_mean_sl_dry, 1),
"TreeGrass" = round(tree_grass_mean_sl_dry, 1),
"Tentative" = round(tentative_mean_sl, 1))
# Construct dynamic labels
labels <- paste(names(mean_values), ": Mean SL = ", mean_values, sep = "")
# Plot
ggplot(sl_dry_mvmt, aes(x = x, y = value, color = name, linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black"),
labels = labels) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted"),
labels = labels) +
theme_bw() +
theme(legend.position = c(0.65, 0.65))Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
Creating von Mises distributions
Code
ta_dry_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))
# Tentative distribution
ta_dry_mvmt$Tentative <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = tentative_kappa)
# Shrub grass
ta_dry_mvmt$ShrubGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = shrub_grass_mvmt_dry$updated_kappa)
# Bare ground
ta_dry_mvmt$Floodplain <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = floodplain_mvmt_dry$updated_kappa)
# Dry grassland
ta_dry_mvmt$DryGrassland <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = dry_grassland_mvmt_dry$updated_kappa)
# Open forest
ta_dry_mvmt$OpenForest <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = open_forest_mvmt_dry$updated_kappa)
# Tree grass
ta_dry_mvmt$TreeGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = tree_grass_mvmt_dry$updated_kappa)
# Pivot from wide to long data
ta_dry_mvmt <- ta_dry_mvmt %>%
pivot_longer(cols = -x)
ta_dry_mvmt$name <- factor(ta_dry_mvmt$name, levels = order)
# Plot
ggplot(ta_dry_mvmt,
aes(x = x, y = value, color = factor(name), linetype = factor(name))
) +
geom_line(alpha = 0.75, linewidth = 1) +
xlab("Turning angle (radians)") +
scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black")) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted")) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Dry season") +
theme_bw() +
theme(legend.position = "bottom")Code
coef_df_wet_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
"scaled_coefficient" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"])
# step length coefficients
sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate))
sl_wet_mvmt_mvparams <- sl_wet_mvmt_mvparams %>% mutate(
scale_factor = sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# log step length coefficients
log_sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate))
log_sl_wet_mvmt_mvparams <- log_sl_wet_mvmt_mvparams %>% mutate(
scale_factor = log_sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# turning angle coefficients
cos_ta_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate))
cos_ta_wet_mvmt_mvparams <- cos_ta_wet_mvmt_mvparams %>% mutate(
scale_factor = cos_ta_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# shrub_grass
shrub_grass_mvmt_wet <- data.frame(
"pheno" = "ShrubGrass",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1],
"updated_scale" = 1/((1/tentative_scale) - sl_wet_mvmt_mvparams$nat_coef[1]),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1]
)
shrub_grass_mean_sl_wet <- shrub_grass_mvmt_wet$updated_shape * shrub_grass_mvmt_wet$updated_scale
# floodplain
floodplain_mvmt_wet <- data.frame(
"pheno" = "Floodplain",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[2],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[2])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[2]
)
floodplain_mean_sl_wet <- floodplain_mvmt_wet$updated_shape * floodplain_mvmt_wet$updated_scale
# dry_grassland
dry_grassland_mvmt_wet <- data.frame(
"pheno" = "wetGrassland",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[3],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[3])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[3]
)
dry_grassland_mean_sl_wet <- dry_grassland_mvmt_wet$updated_shape * dry_grassland_mvmt_wet$updated_scale
# open_forest
open_forest_mvmt_wet <- data.frame(
"pheno" = "OpenForest",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[4],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[4])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[4]
)
open_forest_mean_sl_wet <- open_forest_mvmt_wet$updated_shape * open_forest_mvmt_wet$updated_scale
# tree_grass
tree_grass_mvmt_wet <- data.frame(
"pheno" = "TreeGrass",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[5],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[5])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[5]
)
tree_grass_mean_sl_wet <- tree_grass_mvmt_wet$updated_shape * tree_grass_mvmt_wet$updated_scaleCreating Gamma distributions
Code
sl_wet_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))
# Tentative distribution
sl_wet_mvmt$Tentative <- dgamma(x = sl_wet_mvmt$x,
shape = tentative_shape,
scale = tentative_scale)
# Shrub grass
sl_wet_mvmt$ShrubGrass <- dgamma(x = sl_wet_mvmt$x,
shape = shrub_grass_mvmt_wet$updated_shape,
scale = shrub_grass_mvmt_wet$updated_scale)
# Bare ground
sl_wet_mvmt$Floodplain <- dgamma(x = sl_wet_mvmt$x,
shape = floodplain_mvmt_wet$updated_shape,
scale = floodplain_mvmt_wet$updated_scale)
# Dry grassland
sl_wet_mvmt$DryGrassland <- dgamma(x = sl_wet_mvmt$x,
shape = dry_grassland_mvmt_wet$updated_shape,
scale = dry_grassland_mvmt_wet$updated_scale)
# Open forest
sl_wet_mvmt$OpenForest <- dgamma(x = sl_wet_mvmt$x,
shape = open_forest_mvmt_wet$updated_shape,
scale = open_forest_mvmt_wet$updated_scale)
# Tree grass
sl_wet_mvmt$TreeGrass <- dgamma(x = sl_wet_mvmt$x,
shape = tree_grass_mvmt_wet$updated_shape,
scale = tree_grass_mvmt_wet$updated_scale)
# Pivot from wide to long data
sl_wet_mvmt <- sl_wet_mvmt %>%
pivot_longer(cols = -x)
sl_wet_mvmt$name <- factor(sl_wet_mvmt$name, levels = order)
mean_values_wet <- c("ShrubGrass" = round(shrub_grass_mean_sl_wet, 1),
"Floodplain" = round(floodplain_mean_sl_wet, 1),
"DryGrassland" = round(dry_grassland_mean_sl_wet, 1),
"OpenForest" = round(open_forest_mean_sl_wet, 1),
"TreeGrass" = round(tree_grass_mean_sl_wet, 1),
"Tentative" = round(tentative_mean_sl, 1))
# Construct dynamic labels
labels <- paste(names(mean_values_wet), ": Mean SL = ", mean_values_wet, sep = "")
# Plot
ggplot(sl_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black"),
labels = labels) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted"),
labels = labels) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Wet season") +
theme_bw() +
theme(legend.position = c(0.65, 0.65))Creating von Mises distributions
Code
ta_wet_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))
# Tentative distribution
ta_wet_mvmt$Tentative <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = tentative_kappa)
# Shrub grass
ta_wet_mvmt$ShrubGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = shrub_grass_mvmt_wet$updated_kappa)
# Bare ground
ta_wet_mvmt$Floodplain <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = floodplain_mvmt_wet$updated_kappa)
# Dry grassland
ta_wet_mvmt$DryGrassland <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = dry_grassland_mvmt_wet$updated_kappa)
# Open forest
ta_wet_mvmt$OpenForest <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = open_forest_mvmt_wet$updated_kappa)
# Tree grass
ta_wet_mvmt$TreeGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = tree_grass_mvmt_wet$updated_kappa)
# Pivot from wide to long data
ta_wet_mvmt <- ta_wet_mvmt %>%
pivot_longer(cols = -x)
ta_wet_mvmt$name <- factor(ta_wet_mvmt$name, levels = order)
# Plot
ggplot(ta_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
xlab("Turning angle (radians)") +
scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black")) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted")) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Dry season") +
theme_bw() +
theme(legend.position = "bottom")Session Information
R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Australia/Brisbane
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] amt_0.2.2.0 circular_0.5-1 lemon_0.5.0
[4] beepr_2.0 tictoc_1.2.1 glmmTMB_1.1.11
[7] sjPlot_2.8.17 TwoStepCLogit_1.2.6 lubridate_1.9.4
[10] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[13] purrr_1.0.4 readr_2.1.5 tidyr_1.3.1
[16] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] sjlabelled_1.2.0 tidyselect_1.2.1 farver_2.1.2
[4] fastmap_1.2.0 bayestestR_0.15.3 sjstats_0.19.0
[7] digest_0.6.37 timechange_0.3.0 lifecycle_1.0.4
[10] sf_1.0-20 survival_3.8-3 magrittr_2.0.3
[13] compiler_4.5.0 rlang_1.1.6 tools_4.5.0
[16] yaml_2.3.10 knitr_1.50 labeling_0.4.3
[19] bit_4.6.0 classInt_0.4-11 plyr_1.8.9
[22] RColorBrewer_1.1-3 KernSmooth_2.23-26 withr_3.0.2
[25] numDeriv_2016.8-1.1 grid_4.5.0 datawizard_1.1.0
[28] e1071_1.7-16 scales_1.4.0 MASS_7.3-65
[31] insight_1.2.0 cli_3.6.5 mvtnorm_1.3-3
[34] crayon_1.5.3 rmarkdown_2.29 ragg_1.4.0
[37] reformulas_0.4.1 generics_0.1.3 rstudioapi_0.17.1
[40] performance_0.13.0 tzdb_0.5.0 parameters_0.25.0
[43] minqa_1.2.8 DBI_1.2.3 proxy_0.4-27
[46] audio_0.1-11 splines_4.5.0 parallel_4.5.0
[49] effectsize_1.0.0 vctrs_0.6.5 boot_1.3-31
[52] Matrix_1.7-3 jsonlite_2.0.0 hms_1.1.3
[55] bit64_4.6.0-1 systemfonts_1.2.3 units_0.8-7
[58] glue_1.8.0 nloptr_2.2.1 stringi_1.8.7
[61] gtable_0.3.6 ggeffects_2.2.1 lme4_1.1-37
[64] pillar_1.10.2 htmltools_0.5.8.1 R6_2.6.1
[67] TMB_1.9.17 textshaping_1.0.1 Rdpack_2.6.4
[70] vroom_1.6.5 evaluate_1.0.3 lattice_0.22-6
[73] rbibutils_2.3 class_7.3-23 Rcpp_1.0.14
[76] gridExtra_2.3 nlme_3.1-168 mgcv_1.9-1
[79] xfun_0.52 sjmisc_2.8.10 pkgconfig_2.0.3